/*
 * NewtonSolver.cpp
 *
 *  Created on: 20 Oct 2011
 *      Author: allan
 */

#include "NewtonSolver.h"
#include <Eigen/LU>

NewtonSolver::NewtonSolver() :
tolerance(1.0E-6), maxiterations(20)
{}

NewtonSolver::NewtonSolver(unsigned dimension) :
tolerance(1.0E-6), maxiterations(20)
{
	Initialiaze(dimension);
}

void NewtonSolver::Initialiaze(unsigned dimension)
{
	F = VectorXd::Zero(dimension);
	J = MatrixXd::Zero(dimension, dimension);
}

void NewtonSolver::SetSolverOptions(double tolerance, unsigned maxiterations)
{
	this->tolerance = tolerance;
	this->maxiterations = maxiterations;
}

bool NewtonSolver::Iterate(NewtonProblem& problem, VectorXd& x)
{
	problem.Function(x, F);
	problem.Jacobian(x, J);
	
	FullPivLU<MatrixXd> lu(J);
	
	VectorXd dx = lu.solve(-F);
	
	x += dx;
	
	return F.norm() >= tolerance;
}

bool NewtonSolver::Solve(NewtonProblem& problem, VectorXd& x)
{
	unsigned counter = 0;
	
	while(Iterate(problem, x) && counter++ < maxiterations);
	
	return counter < maxiterations;
}
